<html>
<head>
  <meta HTTP-EQUIV="Content-Type" CONTENT="text/html;charset=ISO-8859-1">
  <title>emgmm.m</title>
<link rel="stylesheet" type="text/css" href="../../../m-syntax.css">
</head>
<body>
<code>
<span class=defun_kw>function</span>&nbsp;<span class=defun_out>model</span>=<span class=defun_name>emgmm</span>(<span class=defun_in>X,options,init_model</span>)<br>
<span class=h1>%&nbsp;EMGMM&nbsp;Expectation-Maximization&nbsp;Algorithm&nbsp;for&nbsp;Gaussian&nbsp;mixture&nbsp;model.</span><br>
<span class=help>%&nbsp;</span><br>
<span class=help>%&nbsp;<span class=help_field>Synopsis:</span></span><br>
<span class=help>%&nbsp;&nbsp;model&nbsp;=&nbsp;emgmm(X)</span><br>
<span class=help>%&nbsp;&nbsp;model&nbsp;=&nbsp;emgmm(X,options)</span><br>
<span class=help>%&nbsp;&nbsp;model&nbsp;=&nbsp;emgmm(X,options,init_model)</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;<span class=help_field>Description:</span></span><br>
<span class=help>%&nbsp;&nbsp;This&nbsp;function&nbsp;implements&nbsp;the&nbsp;Expectation-Maximization&nbsp;algorithm&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;(EM)&nbsp;[Schles68][DLR77]&nbsp;which&nbsp;computes&nbsp;the&nbsp;maximum-likelihood&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;estimate&nbsp;of&nbsp;the&nbsp;paramaters&nbsp;of&nbsp;the&nbsp;Gaussian&nbsp;mixture&nbsp;model&nbsp;(GMM).&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;The&nbsp;EM&nbsp;algorithm&nbsp;is&nbsp;an&nbsp;iterative&nbsp;procedure&nbsp;which&nbsp;monotonically&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;increases&nbsp;log-likelihood&nbsp;of&nbsp;the&nbsp;current&nbsp;estimate&nbsp;until&nbsp;it&nbsp;reaches&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;a&nbsp;local&nbsp;optimum.&nbsp;</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;&nbsp;The&nbsp;number&nbsp;of&nbsp;components&nbsp;of&nbsp;the&nbsp;GMM&nbsp;is&nbsp;given&nbsp;in&nbsp;options.ncomp&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;(default&nbsp;2).</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;&nbsp;The&nbsp;following&nbsp;three&nbsp;stopping&nbsp;are&nbsp;condition&nbsp;used:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;1.&nbsp;Improvement&nbsp;of&nbsp;the&nbsp;log-likelihood&nbsp;is&nbsp;less&nbsp;than&nbsp;given</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;threshold</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;logL(t+1)&nbsp;&nbsp;-&nbsp;logL(t)&nbsp;&lt;&nbsp;options.eps_logL</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;2.&nbsp;Change&nbsp;of&nbsp;the&nbsp;squared&nbsp;differences&nbsp;of&nbsp;a&nbsp;estimated&nbsp;posteriory&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;probabilities&nbsp;is&nbsp;less&nbsp;than&nbsp;given&nbsp;threshold</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;||alpha(t+1)&nbsp;-&nbsp;alpha(t)||^2&nbsp;&lt;&nbsp;options.eps_alpha</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;3.&nbsp;Number&nbsp;of&nbsp;iterations&nbsp;exceeds&nbsp;given&nbsp;threshold.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;t&nbsp;&gt;=&nbsp;options.tmax&nbsp;</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;&nbsp;The&nbsp;type&nbsp;of&nbsp;estimated&nbsp;covariance&nbsp;matrices&nbsp;is&nbsp;optional:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;options.cov_type&nbsp;=&nbsp;'full'&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;full&nbsp;covariance&nbsp;matrix&nbsp;(default)</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;options.cov_type&nbsp;=&nbsp;'diag'&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;diagonal&nbsp;covarinace&nbsp;matrix</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;cov_options.type&nbsp;=&nbsp;'spherical'&nbsp;spherical&nbsp;covariance&nbsp;matrix</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;&nbsp;The&nbsp;initial&nbsp;model&nbsp;(estimate)&nbsp;is&nbsp;selected:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;1.&nbsp;randomly&nbsp;(options.init&nbsp;=&nbsp;'random')&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;2.&nbsp;using&nbsp;C-means&nbsp;(options.init&nbsp;=&nbsp;'cmeans')</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;3.&nbsp;using&nbsp;the&nbsp;user&nbsp;specified&nbsp;init_model.</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;<span class=help_field>Input:</span></span><br>
<span class=help>%&nbsp;&nbsp;X&nbsp;[dim&nbsp;x&nbsp;num_data]&nbsp;Data&nbsp;sample.</span><br>
<span class=help>%&nbsp;&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;options&nbsp;[struct]&nbsp;Control&nbsp;paramaters:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.ncomp&nbsp;[1x1]&nbsp;Number&nbsp;of&nbsp;components&nbsp;of&nbsp;GMM&nbsp;(default&nbsp;2).</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.tmax&nbsp;[1x1]&nbsp;Maximal&nbsp;number&nbsp;of&nbsp;iterations&nbsp;(default&nbsp;inf).</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.eps_logL&nbsp;[1x1]&nbsp;Minimal&nbsp;improvement&nbsp;in&nbsp;log-likelihood&nbsp;(default&nbsp;0).</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.eps_alpha&nbsp;[1x1]&nbsp;Minimal&nbsp;change&nbsp;of&nbsp;Alphas&nbsp;(default&nbsp;0).</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.cov_type&nbsp;[1x1]&nbsp;Type&nbsp;of&nbsp;estimated&nbsp;covarince&nbsp;matrices&nbsp;(see&nbsp;above).</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.init&nbsp;[string]&nbsp;'random'&nbsp;use&nbsp;random&nbsp;initial&nbsp;model&nbsp;(default);</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;'cmeans'&nbsp;use&nbsp;K-means&nbsp;to&nbsp;find&nbsp;initial&nbsp;model.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.verb&nbsp;[1x1]&nbsp;If&nbsp;1&nbsp;then&nbsp;info&nbsp;is&nbsp;displayed&nbsp;(default&nbsp;0).</span><br>
<span class=help>%&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;init_model&nbsp;[struct]&nbsp;Initial&nbsp;model:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Mean&nbsp;[dim&nbsp;x&nbsp;ncomp]&nbsp;Mean&nbsp;vectors.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Cov&nbsp;[dim&nbsp;x&nbsp;dim&nbsp;x&nbsp;ncomp]&nbsp;Covariance&nbsp;matrices.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Priors&nbsp;[1&nbsp;x&nbsp;ncomp]&nbsp;Weights&nbsp;of&nbsp;mixture&nbsp;components.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Alpha&nbsp;[ncomp&nbsp;x&nbsp;num_data]&nbsp;(optional)&nbsp;Distribution&nbsp;of&nbsp;hidden&nbsp;state.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.t&nbsp;[1x1]&nbsp;(optional)&nbsp;Counter&nbsp;of&nbsp;iterations.</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;<span class=help_field>Output:</span></span><br>
<span class=help>%&nbsp;&nbsp;model&nbsp;[struct]&nbsp;Estimated&nbsp;Gaussian&nbsp;mixture&nbsp;model:</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Mean&nbsp;[dim&nbsp;x&nbsp;ncomp]&nbsp;Mean&nbsp;vectors.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Cov&nbsp;[dim&nbsp;x&nbsp;dim&nbsp;x&nbsp;ncomp]&nbsp;Covariance&nbsp;matrices.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.Prior&nbsp;[1&nbsp;x&nbsp;ncomp]&nbsp;Weights&nbsp;of&nbsp;mixture&nbsp;components.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.t&nbsp;[1x1]&nbsp;Number&nbsp;iterations.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.options&nbsp;[struct]&nbsp;Copy&nbsp;of&nbsp;used&nbsp;options.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;.exitflag&nbsp;[int]&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;...&nbsp;maximal&nbsp;number&nbsp;of&nbsp;iterations&nbsp;was&nbsp;exceeded.</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;1&nbsp;or&nbsp;2&nbsp;...&nbsp;EM&nbsp;has&nbsp;converged;&nbsp;indicates&nbsp;which&nbsp;stopping&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;was&nbsp;used&nbsp;(see&nbsp;above).</span><br>
<span class=help>%&nbsp;&nbsp;</span><br>
<span class=help>%&nbsp;<span class=help_field>Example:</span></span><br>
<span class=help>%&nbsp;Note:&nbsp;if&nbsp;EM&nbsp;algorithm&nbsp;does&nbsp;not&nbsp;converge&nbsp;run&nbsp;it&nbsp;again&nbsp;from&nbsp;different</span><br>
<span class=help>%&nbsp;initial&nbsp;model.</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;EM&nbsp;is&nbsp;used&nbsp;to&nbsp;estimate&nbsp;parameters&nbsp;of&nbsp;mixture&nbsp;of&nbsp;2&nbsp;Guassians:</span><br>
<span class=help>%&nbsp;&nbsp;true_model&nbsp;=&nbsp;struct('Mean',[-2&nbsp;2],'Cov',[1&nbsp;0.5],'Prior',[0.4&nbsp;0.6]);</span><br>
<span class=help>%&nbsp;&nbsp;sample&nbsp;=&nbsp;gmmsamp(true_model,&nbsp;100);</span><br>
<span class=help>%&nbsp;&nbsp;estimated_model&nbsp;=&nbsp;emgmm(sample.X,struct('ncomp',2,'verb',1));</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;&nbsp;figure;&nbsp;ppatterns(sample.X);</span><br>
<span class=help>%&nbsp;&nbsp;h1=pgmm(true_model,struct('color','r'));</span><br>
<span class=help>%&nbsp;&nbsp;h2=pgmm(estimated_model,struct('color','b'));</span><br>
<span class=help>%&nbsp;&nbsp;legend([h1(1)&nbsp;h2(1)],'Ground&nbsp;truth',&nbsp;'ML&nbsp;estimation');&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;figure;&nbsp;hold&nbsp;on;&nbsp;xlabel('iterations');&nbsp;ylabel('log-likelihood');</span><br>
<span class=help>%&nbsp;&nbsp;plot(&nbsp;estimated_model.logL&nbsp;);</span><br>
<span class=help>%</span><br>
<span class=help>%&nbsp;See&nbsp;also&nbsp;</span><br>
<span class=help>%&nbsp;&nbsp;MLCGMM,&nbsp;MMGAUSS,&nbsp;PDFGMM,&nbsp;GMMSAMP.</span><br>
<span class=help>%</span><br>
<hr>
<span class=help1>%&nbsp;<span class=help1_field>About:</span>&nbsp;Statistical&nbsp;Patte7rn&nbsp;Recognition&nbsp;Toolbox</span><br>
<span class=help1>%&nbsp;(C)&nbsp;1999-2003,&nbsp;Written&nbsp;by&nbsp;Vojtech&nbsp;Franc&nbsp;and&nbsp;Vaclav&nbsp;Hlavac</span><br>
<span class=help1>%&nbsp;&lt;a&nbsp;href="http://www.cvut.cz"&gt;Czech&nbsp;Technical&nbsp;University&nbsp;Prague&lt;/a&gt;</span><br>
<span class=help1>%&nbsp;&lt;a&nbsp;href="http://www.feld.cvut.cz"&gt;Faculty&nbsp;of&nbsp;Electrical&nbsp;Engineering&lt;/a&gt;</span><br>
<span class=help1>%&nbsp;&lt;a&nbsp;href="http://cmp.felk.cvut.cz"&gt;Center&nbsp;for&nbsp;Machine&nbsp;Perception&lt;/a&gt;</span><br>
<br>
<span class=help1>%&nbsp;<span class=help1_field>Modifications:</span></span><br>
<span class=help1>%&nbsp;26-jul-07,&nbsp;VF,&nbsp;inconsistent&nbsp;parameter&nbsp;names&nbsp;'ncomp&nbsp;and&nbsp;'num_gauss'&nbsp;removed</span><br>
<span class=help1>%&nbsp;26-may-2004,&nbsp;VF,&nbsp;initialization&nbsp;by&nbsp;K-means&nbsp;added</span><br>
<span class=help1>%&nbsp;1-may-2004,&nbsp;VF</span><br>
<span class=help1>%&nbsp;19-sep-2003,&nbsp;VF</span><br>
<span class=help1>%&nbsp;16-mar-2003,&nbsp;VF</span><br>
<br>
<br>
<hr>
<span class=comment>%&nbsp;processing&nbsp;input&nbsp;arguments&nbsp;</span><br>
<span class=comment>%&nbsp;-----------------------------------------</span><br>
<span class=keyword>if</span>&nbsp;<span class=stack>nargin</span>&nbsp;&lt;&nbsp;2,&nbsp;options=[];&nbsp;<span class=keyword>else</span>&nbsp;options=c2s(options);&nbsp;<span class=keyword>end</span><br>
<br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'ncomp'</span>),&nbsp;options.ncomp&nbsp;=&nbsp;2;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'tmax'</span>),&nbsp;options.tmax&nbsp;=&nbsp;inf;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'eps_alpha'</span>),&nbsp;options.eps_alpha&nbsp;=&nbsp;0;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'eps_logL'</span>),&nbsp;options.eps_logL&nbsp;=&nbsp;0;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'cov_type'</span>),&nbsp;options.cov_type&nbsp;=&nbsp;<span class=quotes>'full'</span>;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'init'</span>),&nbsp;options.init&nbsp;=&nbsp;<span class=quotes>'random'</span>;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>if</span>&nbsp;~isfield(&nbsp;options,&nbsp;<span class=quotes>'verb'</span>),&nbsp;options.verb&nbsp;=&nbsp;0;&nbsp;<span class=keyword>end</span><br>
<br>
[dim,num_data]&nbsp;=&nbsp;size(X);<br>
<br>
<span class=comment>%&nbsp;setup&nbsp;initial&nbsp;model&nbsp;</span><br>
<span class=comment>%&nbsp;---------------------------------</span><br>
<br>
<span class=keyword>if</span>&nbsp;<span class=stack>nargin</span>&nbsp;==&nbsp;3,<br>
<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;take&nbsp;model&nbsp;from&nbsp;input</span><br>
&nbsp;&nbsp;<span class=comment>%-----------------------------</span><br>
<br>
&nbsp;&nbsp;model&nbsp;=&nbsp;init_model;&nbsp;<br>
&nbsp;&nbsp;<span class=keyword>if</span>&nbsp;~isfield(model,<span class=quotes>'t'</span>),&nbsp;model.t&nbsp;=&nbsp;0;&nbsp;<span class=keyword>end</span><br>
&nbsp;&nbsp;<span class=keyword>if</span>&nbsp;~isfield(model,<span class=quotes>'Alpha'</span>),&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;model.Alpha=-inf*ones(options.ncomp,num_data);<br>
&nbsp;&nbsp;<span class=keyword>end</span><br>
&nbsp;&nbsp;<span class=keyword>if</span>&nbsp;~isfield(model,<span class=quotes>'logL'</span>),&nbsp;model.logL=-inf;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>else</span><br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;compute&nbsp;initial&nbsp;model</span><br>
&nbsp;&nbsp;<span class=comment>%------------------------------------</span><br>
&nbsp;&nbsp;<span class=keyword>switch</span>&nbsp;options.init,<br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;random&nbsp;model</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=label>case</span>&nbsp;<span class=quotes>'random'</span>&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;takes&nbsp;randomly&nbsp;first&nbsp;ncomp&nbsp;trn.&nbsp;vectors&nbsp;as&nbsp;mean&nbsp;vectors</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;inx&nbsp;=&nbsp;randperm(num_data);&nbsp;&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;inx=inx(1:options.ncomp);<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;centers_X&nbsp;=&nbsp;X(:,inx);<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;K-means&nbsp;clustering</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=label>case</span>&nbsp;<span class=quotes>'cmeans'</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;tmp&nbsp;=&nbsp;cmeans(&nbsp;X,&nbsp;options.ncomp&nbsp;);<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;centers_X&nbsp;=&nbsp;tmp.X;<br>
<br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=label>otherwise</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<span class=error>error</span>(<span class=quotes>'Unknown&nbsp;initialization&nbsp;method.'</span>);<br>
&nbsp;&nbsp;<span class=keyword>end</span><br>
<br>
&nbsp;&nbsp;knn&nbsp;=&nbsp;knnrule({<span class=quotes>'X'</span>,centers_X,<span class=quotes>'y'</span>,[1:options.ncomp]},1);<br>
&nbsp;&nbsp;y&nbsp;=&nbsp;knnclass(X,knn);<br>
<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;uses&nbsp;ML&nbsp;estimation&nbsp;of&nbsp;complete&nbsp;data</span><br>
&nbsp;&nbsp;model&nbsp;=&nbsp;mlcgmm(&nbsp;{<span class=quotes>'X'</span>,X,<span class=quotes>'y'</span>,y},&nbsp;options.cov_type&nbsp;);<br>
<br>
&nbsp;&nbsp;model.Alpha&nbsp;=&nbsp;zeros(options.ncomp,num_data);<br>
&nbsp;&nbsp;<span class=keyword>for</span>&nbsp;i&nbsp;=&nbsp;1:options.ncomp,<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.Alpha(i,find(y==i))&nbsp;=&nbsp;1;<br>
&nbsp;&nbsp;<span class=keyword>end</span><br>
&nbsp;&nbsp;model.logL=&nbsp;-inf;<br>
&nbsp;&nbsp;model.t&nbsp;=&nbsp;1;<br>
&nbsp;&nbsp;model.options&nbsp;=&nbsp;options;<br>
&nbsp;&nbsp;model.fun&nbsp;=&nbsp;<span class=quotes>'pdfgmm'</span>;<br>
&nbsp;&nbsp;<br>
<span class=keyword>end</span><br>
<br>
<br>
<span class=comment>%&nbsp;Main&nbsp;loop&nbsp;of&nbsp;EM&nbsp;algorithm&nbsp;</span><br>
<span class=comment>%&nbsp;-------------------------------------</span><br>
model.exitflag&nbsp;=&nbsp;0;<br>
<span class=keyword>while</span>&nbsp;model.exitflag&nbsp;==&nbsp;0&nbsp;&&nbsp;model.t&nbsp;&lt;&nbsp;options.tmax,<br>
<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;counter&nbsp;of&nbsp;iterations</span><br>
&nbsp;&nbsp;model.t&nbsp;=&nbsp;model.t&nbsp;+&nbsp;1;<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=comment>%----------------------------------------------------</span><br>
&nbsp;&nbsp;<span class=comment>%&nbsp;E-Step</span><br>
&nbsp;&nbsp;<span class=comment>%&nbsp;The&nbsp;distribution&nbsp;of&nbsp;hidden&nbsp;states&nbsp;is&nbsp;computed&nbsp;based</span><br>
&nbsp;&nbsp;<span class=comment>%&nbsp;on&nbsp;the&nbsp;current&nbsp;estimate.</span><br>
&nbsp;&nbsp;<span class=comment>%----------------------------------------------------</span><br>
<br>
&nbsp;&nbsp;newAlpha&nbsp;=&nbsp;(model.Prior(:)*ones(1,num_data)).*pdfgauss(X,&nbsp;model);<br>
&nbsp;&nbsp;newLogL&nbsp;=&nbsp;sum(log(sum(newAlpha,1)));&nbsp;&nbsp;<br>
&nbsp;&nbsp;newAlpha&nbsp;=&nbsp;newAlpha./(ones(options.ncomp,1)*sum(newAlpha,1));<br>
<br>
&nbsp;&nbsp;<span class=comment>%------------------------------------------------------</span><br>
&nbsp;&nbsp;<span class=comment>%&nbsp;Stopping&nbsp;conditions.</span><br>
&nbsp;&nbsp;<span class=comment>%------------------------------------------------------</span><br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;1)&nbsp;change&nbsp;in&nbsp;distribution&nbsp;of&nbsp;hidden&nbsp;state&nbsp;Alpha</span><br>
&nbsp;&nbsp;model.delta_alpha&nbsp;=&nbsp;sum(sum((model.Alpha&nbsp;-&nbsp;newAlpha).^2));<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=comment>%&nbsp;2)&nbsp;change&nbsp;in&nbsp;log-Likelihood</span><br>
&nbsp;&nbsp;model.delta_logL&nbsp;=&nbsp;newLogL&nbsp;-&nbsp;model.logL(<span class=keyword>end</span>);<br>
&nbsp;&nbsp;model.logL&nbsp;=&nbsp;[model.logL&nbsp;newLogL];<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=keyword>if</span>&nbsp;options.verb,<br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=io>fprintf</span>(<span class=quotes>'%d:&nbsp;logL=%f,&nbsp;delta_logL=%f,&nbsp;delta_alpha=%f\n'</span>,...<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;model.t,&nbsp;model.logL(<span class=keyword>end</span>),&nbsp;model.delta_logL,&nbsp;model.delta_alpha&nbsp;);<br>
&nbsp;&nbsp;<span class=keyword>end</span><br>
<br>
&nbsp;&nbsp;<span class=keyword>if</span>&nbsp;options.eps_logL&nbsp;&gt;=&nbsp;model.delta_logL,<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.exitflag&nbsp;=&nbsp;1;<br>
&nbsp;&nbsp;<span class=keyword>elseif</span>&nbsp;options.eps_alpha&nbsp;&gt;=&nbsp;model.delta_alpha,<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.exitflag&nbsp;=&nbsp;2;<br>
&nbsp;&nbsp;<span class=keyword>else</span><br>
<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.Alpha&nbsp;=&nbsp;newAlpha;<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%----------------------------------------------------</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;M-Step</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;The&nbsp;new&nbsp;parameters&nbsp;maximizing&nbsp;expectation&nbsp;of&nbsp;</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%&nbsp;log-likelihood&nbsp;are&nbsp;computed.</span><br>
&nbsp;&nbsp;&nbsp;&nbsp;<span class=comment>%----------------------------------------------------</span><br>
<br>
&nbsp;&nbsp;&nbsp;&nbsp;tmp_model&nbsp;=&nbsp;melgmm(X,model.Alpha,options.cov_type);<br>
&nbsp;&nbsp;&nbsp;&nbsp;<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.Mean&nbsp;=&nbsp;tmp_model.Mean;<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.Cov&nbsp;=&nbsp;tmp_model.Cov;<br>
&nbsp;&nbsp;&nbsp;&nbsp;model.Prior&nbsp;=&nbsp;tmp_model.Prior;<br>
&nbsp;&nbsp;<br>
&nbsp;&nbsp;<span class=keyword>end</span><br>
<span class=keyword>end</span>&nbsp;<span class=comment>%&nbsp;while&nbsp;main&nbsp;loop</span><br>
<br>
<span class=jump>return</span>;<br>
<br>
<br>
</code>
